function b_idat = fsb_get_extraction_mask(imean,imean_avg,threshold,coord)

% FSB : Extract a brain, remove muscle and noise, reduce image size
% without loss of information
% use mean idat and a variable threshold.
%
% EXAMPLE:
% b_idat = fsb_get_extraction_mask(imean,imean_avg,threshold,coord)
%
% INPUT:
% imean:            3D Image data (mean of 4D volume)
% imean_avg:        Average image intensity value
% threshold:        Threshold can be set to values between 0 and about
%                   Values around 1000 work best.
% coord:            struct with coordinates of seed point for extraction
%                   contains subfields slice_x, slice_y, slice_z with
%                   the coordinates
%
% OUTPUT:
% b_idat: An extraction mask for a brain
%
% CALLED BY:
% FSB.m
%
% NOTES:
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% check if seed point is provided
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if nargin==4
    seed = [coord.slice_x coord.slice_y coord.slice_z];
    imean_avg = avg(imean(round(coord.slice_x-coord.slice_x/6):round(coord.slice_x+coord.slice_x/6),...
        round(coord.slice_y-coord.slice_y/6):round(coord.slice_y+coord.slice_y/6),...
        round(coord.slice_z-coord.slice_z/6):round(coord.slice_z+coord.slice_z/6))); % calculate image average value to use as threshold
else
    [a,b,c]=size(imean);
    seed = round([a,b,c]/2);
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check if threshold is set correctly and create threshold image
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



if threshold>imean_avg;
    threshold = imean_avg;
    disp(sprintf(strcat('Threshold too high, reset to max value ( ', num2str(threshold),')\nThis results in an inclusion of all voxels in the image')));
end

    
imask = imean>(imean_avg-threshold);

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% check that seed point is within the brain
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
i = 0;
while imask(seed(1),seed(2),seed(3))<1;
    seed= seed+[1,1,1];
    i = i+1;
    if i == 10;
        disp('please reposition seed point, cannot find brain');
        b_idat = imean;
        return
    end
end

     

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%use external function to fill the brain from inside out
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
B = FloodFill3D_s(imask, seed); % FloodFill3D is copyright Faheem Dinath

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%smooth filled brain to increase size
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ismooth = single(smooth3(B,'box',3));

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%clear variables to save memory
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear B;

ithres = ismooth==0;%use this to generate new threshold image
imean(ithres)=0;
b_idat =imean;

end
